Los modelos jerárquicos involucran varios parámetros de tal manera que las creencias de unos de los parámetros dependen de manera significativa de los valores de otros parámetros de manera jerárquica.

Para ilustrar esta dependencia consideraremos un ejemplo. Consideremos el caso en el que tenemos varias monedas acuñadas en la misma casa de monedas, es razonable pensar que una fábrica sesgada a águilas tenderá a producir monedas con sesgo hacia águilas. La estimación del sesgo de una moneda depende de la estimación del sesgo de la fábrica que a su vez está influido por los datos de todas las monedas. Veremos que la estructura de dependenica a lo largo de los parámetros generan estimaciones mejor informadas para todos los parámetros.

Si pensamos únicamente en dos monedas que provienen de la misma casa de moneda tenemos:

  1. Creencias iniciales de los posibles valores de los parámetros (sesgos de las monedas).

  2. Creencias iniciales de la dependencia de los parámetros por provenir de la misma fábrica.

Cuando observamos lanzamientos de las monedas actualizamos nuestras creencias relativas a los sesgos de las monedas y también actualizamos nuestras creencias acerca de la dependencia de los parámetros.

Ejemplo: Modelo jerárquico de una moneda

Recordemos el caso de lanzamientos de una moneda. Le asignamos una distribución inicial beta, recordemos que la distribución beta tienen dos parámetors \(a\) y \(b\):

\[p(\theta)=\frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1}\]

Con el fin de hacer los parámetros más intuitivos los podemos expresar en términos de la media \(\mu\) de la distribución inicial y el tamaño de muestra \(K\), el cual refleja la variabilidad de la distribución, esto es, la confianza que tenemos en la media. Los parámetros correspondientes en la distribución beta son:

\[a=\mu K, \quad b = (1-\mu)K\]

Ahora introducimos el modelo jerárquico. En lugar de especificar un valor particular para \(\mu\), consideramos que \(\mu\) puede tomar distintos valores (entre 0 y 1), y definimos una distribución de probabilidad sobre esos valores. Podemos pensar que esta distribución describe nuestra incertidumbre acerca de la construcción de la máquina que manufactura las monedas.

Veremos posteriormente que en el caso de más de una moneda el modelo permite que cada moneda tenga un sesgo distinto, pero ambas tenderán a tener un sesgo cercano a \(\mu\), algunas aleatoriamente tendrán un valor de \(\theta\) mayor a \(\mu\) y otras menor. Entre más grande \(K\) mayor será la consistencia de la acuñadora y los valores \(\theta\) serán más cercanos a \(\mu\). Si observamos varios lanzamientos de una moneda tendremos información tanto de \(\theta\) como de \(\mu\).

Para hacer un análisis bayesiano aún nos hace falta definir la distribución inicial sobre el parámetro \(\mu\); utilicemos una distribución Beta:

\[p(\mu)=Beta(\mu|A_{\mu}, B_{\mu})\]

donde \(A_{\mu}\) y \(B_{\mu}\) se conocen como hiperparámetros y son constantes. En este caso estamos consideramos que \(\mu\) se ubica típicamente cerca de la media \(A_{\mu}/(A_{\mu} + B_{\mu})\) y tomaremos \(K\) como constante.

Recordemos también que en el ejemplo de una moneda teníamos que la función de verosimilitud es una distribución Bernoulli para los datos:

\[p(x|\theta) = \theta^x(1-\theta)^{1-x}\]

Si utilizamos las iniciales Beta para \(\mu\) y \(\theta\) como discutimos arriba, solo nos resta aplicar la regla de Bayes con nuestros dos parámetros desconocidos \(\mu\) y \(\theta\):

\[p(\theta, \mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\]

Hay dos aspectos a considerar en el problema:

  1. La verosimilitud no depende de \(\mu\) por lo que \[p(x|\theta, \mu)=p(x|\theta)\]

  2. La distribución inicial en el espacio de parámetros bivariado \((\theta,\mu)\) se puede factorizar como sigue:

\[p(\theta,\mu)=p(\theta|\mu)p(\mu)\]

Por lo tanto \[p(\theta,\mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\] \[=\frac{p(x|\theta)p(\theta|\mu)p(\mu)}{p(x)}\] El siguiente modelo gráfico resume las independencias condicionales de la última ecuación:

library(DiagrammeR)
#> Warning: package 'DiagrammeR' was built under R version 3.5.1
grViz("
digraph betabinom {

  # a 'graph' statement
  graph [overlap = true]

  node [shape = box,penwidth = 0.1, fixedsize = true, fontsize = 9,
        fontname = Helvetica, width = 0.2, height = 0.2]
  AB[label = 'A,B']; 

  node [shape = circle,
        fixedsize = true, fontsize = 9,
        fontname = Helvetica, width = 0.2] // sets as circles
  mu[label = 'μ']; theta[label = 'θ']; x;

  # several 'edge' statements
  edge[color = grey, arrowsize = 0.5, penwidth = 0.5]
  AB->mu; mu->theta; theta->x;
}
")

Aproximación por cuadrícula

En el caso jerárquico, no se puede derivar la distribución posterior de manera analítica pero si los parámetros e hiperparámetros toman un número finito de
valores y no hay muchos de ellos, podemos aproximar la posterior usando aproximación por cuadrícula.

A continuación graficamos las distribuciones correspondientes al caso en que la distribución inicial del hiperparámetro \(\mu\) tiene la forma de una distribución \(Beta(2, 2)\), es decir, creemos que la media de la acuñadora \(\mu\) es 0.5, pero existe bastante incertidumbre acerca del valor central.

La distribución de \(\theta\), esto es, la distribución inicial que refleja la dependencia entre \(\theta\) y \(\mu\), se expresa por medio de otra distribución beta; en el ejemplo utilizaremos el valor fijo K=100:

\[\theta|\mu \sim Beta(\mu 100, (1-\mu)100)\]

Esta distribución inicial expresa un alto grado de certeza que una acuñadora con hiperparámetro \(\mu\) genera monedas con sesgo cercano a \(\mu\):

A_mu <- 2; B_mu <- 2; K <- 100

# Distribución inicial conjunta p(theta, mu)
p_conjunta <- function(mu, theta, A_mu, B_mu, K){
  # marginal p(mu)
  p_mu <- dbeta(mu, A_mu, B_mu)
  # condicional p(theta | mu)
  p_theta_mu <- dbeta(theta, mu * K, (1 - mu) * K)
  p_mu * p_theta_mu
}

grid <- expand.grid(theta = seq(0.01, 0.99, 0.005), 
  mu = seq(0.01, 0.99, 0.005))

grid_inicial <- grid %>% 
  mutate(p_inicial = p_conjunta(mu, theta, A_mu, B_mu, K))

plot_conj <- ggplot(grid_inicial, aes(x = theta, y = mu, z = p_inicial)) + 
  stat_contour(binwidth = 2, aes(color = ..level..)) +
  scale_x_continuous(expression(theta), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta,mu)), guide = FALSE) 

grid_mu <- grid_inicial %>%
  group_by(mu) %>%
  summarise(p_mu = sum(p_inicial)) %>%
  ungroup() %>%
  mutate(p_mu = p_mu / sum(p_mu))

plot_mu <- ggplot(grid_mu, aes(x = mu, y = p_mu)) +
  geom_path() +
  labs(y = expression(p(mu)), 
       x = expression(mu)) + 
  ylim(0, 0.015) + coord_flip()

mus <- rbeta(5000, A_mu, B_mu)
sims_marg <- data.frame(sims = 1:5000, 
  sims_marg = rbeta(5000, mus * 100, (1 - mus) * 100))

plot_theta <- ggplot(sims_marg, aes(x = sims_marg)) + 
  geom_density(adjust = 2) +
  labs(y = expression(p(theta)), 
       x = expression(theta))

# p(theta|mu=0.75)
sims_cond_1 <- rbeta(5000, 0.75 * 100, (1 - 0.75) * 100)
# p(theta|mu=0.25)
sims_cond_2 <- rbeta(5000, 0.25 * 100, (1 - 0.25) * 100)

# marginal = sims_marg

sims <- data.frame(sim = 1:5000, 
  cond_1 = sims_cond_1, 
  cond_2 = sims_cond_2) %>%
  gather(dist, values, -sim) %>%
  mutate(mu = ifelse(dist == "cond_1", 0.75, 0.25))


plots_marginales <- ggplot(sims, aes(x = values)) +
  labs(y = expression(p(paste(theta, l, mu))), 
       x = expression(theta)) + 
  geom_density(adjust = 2) + 
  facet_wrap(~ mu, ncol = 1)

grid.arrange(plot_conj, plot_mu, plot_theta, plots_marginales, ncol=2)

Verosimilitud.


likeBern <- function(z, N){
  function(theta){
    theta ^ z * (1 - theta) ^ (N - z)
  }
}

# Valores observados
z <- 9; N <- 12
mi_like <- likeBern(z, N)

grid_like <- expand.grid(x = seq(0.01, 1, 0.05), y = seq(0.01, 1, 0.05))
grid_like <- grid_like %>% 
  mutate(z = mi_like(x))

ggplot(grid_like, aes(x = x, y = y, z = z)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(L(theta,mu)), guide = FALSE) +
  coord_fixed()


# cuadrícula
grid_post <- grid_inicial %>%
  mutate(
    prior = p_inicial / sum(p_inicial),
    Like = theta ^ z * (1 - theta) ^ (N - z), # verosimilitud 
    posterior = (Like * prior) / sum(Like * prior)
  )  

post_conj <- ggplot(grid_post, aes(x = theta, y = mu, z = posterior)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta,mu)), guide = FALSE)

Posterior.

# head(grid_post)

grid_mu <- grid_post %>%
  group_by(mu) %>%
  summarise(post_mu = sum(posterior)) %>%
  ungroup() %>%
  mutate(post_mu = post_mu / sum(post_mu))

post_mu <- ggplot(grid_mu, aes(x = mu, y = post_mu)) +
  geom_path() +
  labs(y = expression(p(paste(mu, l, x))), 
       x = expression(mu))+ coord_flip()

grid_theta <- grid_post %>%
  group_by(theta) %>%
  summarise(post_theta = sum(posterior)) %>%
  ungroup() %>%
  mutate(post_theta = post_theta / sum(post_theta))

post_theta <- ggplot(grid_theta, aes(x = theta, y = post_theta)) +
  geom_path() +
  labs(y = expression(p(paste(theta, l, x))), 
       x = expression(theta))

grid_theta_m <- grid_post %>%
  filter(mu == 0.75 | mu == 0.25) %>%
  group_by(theta, mu) %>%
  summarise(post_theta = sum(posterior)) %>%
  group_by(mu) %>%
  mutate(post_theta = post_theta / sum(post_theta))

post_marginales <- ggplot(grid_theta_m, aes(x = theta, y = post_theta)) +
  geom_path() + 
  facet_wrap(~ mu, ncol = 1) +
  labs(y = expression(p(paste(theta, l, x, mu))), 
       x = expression(theta))

Ejemplo: Múltiples monedas de una misma fábrica

La sección anterior considera el escenario en que lanzamos una moneda y hacemos inferencia del parámetro de sesgo \(\theta\) y del hiperparámetro \(\mu\). Ahora consideramos recolectar información de múltiples monedas, si cada moneda tiene su propio sesgo \(\theta_j\) entonces debemos estimar un parámetro distinto para cada moneda.

Suponemos que todas las monedas provienen de la misma fábrica, esto implica que tenemos la mismas creencias iniciales \(\mu\) para todas las monedas. Suponemos también que cada moneda se acuña de manera independiente, esto es, que condicional al parámetro \(\mu\) los parámetros \(\theta_j\) son independientes en nuestras creencias iniciales.

Aproximación por cuadrícula

Supongamos que tenemos dos monedas de la misma fábrica. El objetivo es estimar los sesgos \(\theta_1\), \(\theta_2\) de las dos monedas y estimar simultáneamente el parámetro \(\mu\) correspondiente a la casa de moneda que las fabricó.

Inicial. Asumimos la misma distribución inicial para \(\theta_1\) y \(\theta_2\).

Verosimilitud.


z_1 <- 3; N_1 <- 15
z_2 <- 4; N_2 <- 5

mi_like_1 <- likeBern(z_1, N_1)
mi_like_2 <- likeBern(z_2, N_2)

grid_like <- grid_like %>% 
  mutate(L_1 = mi_like_1(x), 
         L_2 = mi_like_2(x))

plot_like_1 <- ggplot(grid_like, aes(x = x, y = y, z = L_1)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(L(theta,mu)), guide = FALSE)

plot_like_2 <- ggplot(grid_like, aes(x = x, y = y, z = L_2)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(L(theta,mu)), guide = FALSE)

grid.arrange(plot_like_1, plot_like_2, ncol = 3)

Posterior.

JAGS

La sección anterior utilizó un modelo simplificado con el objetivo de poder visualizar el espacio de parámetros. Ahora incluiremos más parámetros con el objetivo de hacer el problema más realista. En los ejemplos anteriores fijamos el grado de dependencia entre \(\mu\) y \(\theta\) de manera arbitraria, a través de \(K\), de tal manera que si \(K\) era grande, los valores \(\theta_j\) individuales se situaban cerca de \(\mu\) y cuando \(K\) era pequeña se permitía más variación.

En situaciones reales no conocemos el valor de \(K\) por adelantado, por lo que dejamos que los datos nos informen acerca de valores creíbles para \(K\). Intuitivamente, cuando la proporción de águilas observadas es similar a lo largo de las monedas, tenemos evidencia de que \(K\) es grande, mientras que si estas proporciones difieren mucho, tenemos evidencia de que \(K\) es pequeña. Debido a que \(K\) pasará de ser una constante a ser un parámetro lo llamaremos \(\kappa\).

La distribución inicial para \(\kappa\) será una Gamma. Elegimos \(\kappa \sim Gamma(1, 0.1)\); ésta tiene media 10 y desviación estándar 10.

library(R2jags)
#> Warning: package 'R2jags' was built under R version 3.5.1
#> Loading required package: rjags
#> Loading required package: coda
#> Linked to JAGS 4.3.0
#> Loaded modules: basemod,bugs
#> 
#> Attaching package: 'R2jags'
#> The following object is masked from 'package:coda':
#> 
#>     traceplot

modelo_jer.txt <-
'
model{
  for(t in 1:N){
    x[t] ~ dbern(theta[coin[t]])
  }
  for(j in 1:nCoins){
    theta[j] ~ dbeta(a, b)
  }
  a <- mu * kappa
  b <- (1 - mu) * kappa
  # A_mu = 2, B_mu = 2
  mu ~ dbeta(2, 2)
  kappa ~ dgamma(1, 0.1)
}
'

Los datos consisten en 5 monedas, cada una se lanza 5 veces cada una, resultando 4 de ellas en 1 águila y 4 soles y otra en 3 águilas y 2 soles.

cat(modelo_jer.txt, file = 'modelo_jer.bugs')

x <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1)
coin <- c(rep(1, 5), rep(2, 5), rep(3, 5), rep(4, 5), rep(5, 5))

jags.inits <- function(){
  list("mu" = runif(1, 0.1, 0.9),
    "kappa" = runif(1, 5, 20))
}

jags_fit <- jags(
  model.file = "modelo_jer.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = x, coin = coin, nCoins = 5,  N = 25),    # lista con los datos
  parameters.to.save = c("mu", "kappa", "theta"),  # parámetros por guardar
  n.chains = 3,   # número de cadenas
  n.iter = 10000,    # número de pasos
  n.burnin = 1000   # calentamiento de la cadena
  )
#> module glm loaded
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 25
#>    Unobserved stochastic nodes: 7
#>    Total graph size: 65
#> 
#> Initializing model

traceplot(jags_fit, varname = c("kappa", "mu", "theta"))


jags_fit
#> Inference for Bugs model at "modelo_jer.bugs", fit using jags,
#>  3 chains, each with 10000 iterations (first 1000 discarded), n.thin = 9
#>  n.sims = 3000 iterations saved
#>          mu.vect sd.vect   2.5%    25%    50%    75%  97.5% Rhat n.eff
#> kappa     14.874  10.838  2.213  6.824 12.164 19.819 42.307    1  1500
#> mu         0.335   0.098  0.163  0.267  0.328  0.400  0.537    1  3000
#> theta[1]   0.287   0.128  0.070  0.195  0.277  0.369  0.560    1  3000
#> theta[2]   0.291   0.129  0.067  0.195  0.283  0.377  0.559    1  3000
#> theta[3]   0.290   0.130  0.069  0.196  0.281  0.373  0.573    1  3000
#> theta[4]   0.292   0.127  0.073  0.200  0.281  0.377  0.562    1  3000
#> theta[5]   0.419   0.144  0.165  0.319  0.408  0.512  0.733    1  3000
#> deviance  30.354   2.124 27.481 28.873 29.940 31.399 35.591    1  3000
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 2.3 and DIC = 32.6
#> DIC is an estimate of expected predictive error (lower deviance is better).

Podemos hacer histogramas de las distribuciones posteriores de los parámetros:

sims_df <- data.frame(n_sim = 1:jags_fit$BUGSoutput$n.sims,
  jags_fit$BUGSoutput$sims.matrix) %>% 
  dplyr::select(-deviance) %>%
  gather(parametro, value, -n_sim)

ggplot(filter(sims_df, parametro != "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8) + facet_wrap(~ parametro)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.



ggplot(filter(sims_df, parametro == "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Referencias